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CONFIDENCE  INTERVAL  ESTIMATION  FOR  SIMULATION  RESPONSE  SURFACES 


by 

Michael  A.  Crane 


1.  Introduction  and  Summary 

In  this  paper,  we  continue  the  investigation  of  response- surface  methodology- 
in  simulation,  begun  in  [2],  [5],  and  [6].  In  [2]  we  showed  that  if  g is  a 
linear  output  function  of  an  input  parameter  A,  over  some  range  a < A < b, 
and  if  statistical  confidence  intervals  for  g(A)  are  obtained  at  measurement 
points  A = A^  and  A = A^,  with  a < A^  < A^  < b,  then  a confidence  band  can 
be  obtained  for  the  entire  function  g over  the  interval  a < A < b without 
making  additional  measurements.  In  [3],  this  result  was  extended  for  the  case 
where  g is  a quadratic  function  of  one  variable,  and  in  [6] , g was  assumed 
to  be  a general  order  polynomial  of  several  variables. 

One  difficulty  arising  in  the  above  methodology  is  that  the  width  of  the 
confidence  band  obtained  is  somewhat  wider  than  would  be  necessary  if  one 
only  required  confidence  intervals  at  individual  points  rather  than  a confidence 
band  for  the  entire  function.  We  are  thus  motivated  to  consider  the  problem  of 
interval  estimation,  as  opposed  to  band  estimation  and  we  take  up  this  subject 
in  this  paper.  We  do  so  for  the  case  of  a linear  function  of  one  variable. 
Generalizations  to  other  types  of  functions  will  be  the  subject  of  subsequent 
papers . 

The  remainder  of  this  paper  is  organized  as  follows.  The  basic  result  is 
derived  in  Section  2.  It  is  shown  that  if  observations  are  taken  at  the  points 
A^  and  A^  and  if  certain  other  conditions  are  satisfied,  then  a confidence 
interval  may  be  obtained  for  g(A)  at  any  point  a < A < b.  In  Section  J>, 


1 


2. 


it  is  shown  that  the  regenerative  method  for  analyzing  simulations  provides  a 
basis  for  applying  the  result  of  Section  2.  An  illustration  is  given  of  a 
simulation  of  an  M/M/l  queue  and  the  results  obtained  are  compared  to  the 
results  obtained  using  the  confidence  band  methodology  of  [2],  [3],  and  [6]. 


2 . The  Result 

Suppose  that  A is  an  input  parameter  for  a simulation  and  that  we  wish 
to  estimate  some  output  parameter  g(A)  which  is  a linear  function  of  A 
over  a range  a < A < b.  Independent  simulation  runs  resulting  in  statistical 
observations  are  made  at  two  particular  parameter  settings  A = A^  and  A = Ag 
where  a < A^  < Ag  < b.  (These  observations  might  be,  for  example,  sample 
variates  from  populations  with  means  gCA^)  and  g(Ag).  Alternatively,  they 
might  be  observations  based  on  random  tours  in  an  application  of  the  regener- 
ative approach;  see  Section  3 for  details  in  this  latter  case.)  At  the 
parameter  setting  A = A^,  the  simulation  produces  n^  statistical  observations 
resulting  in  sample  statistics  r(A^,n^)  and  q(A^,n^).  Similarly,  at  the 
parameter  setting  A = Ag,  ng  observations  result  in  statistics  r(Ag,ng) 
and  ^(A^ng). 

The  following  proposition  forms  a basis  for  a method  to  compute  confidence 
intervals  for  g(A),  where  a < A < b. 

PROPOSITION  1.  Suppose  that  T)(A^,n^)  — > tj(A^)  in  probability  as  n^  — ■>  » 
and  that  q(Ag,ng)  — •>  q(Ag)  in  probability  as  ng  — > ®.  Suppose  further 


that 


3- 


"»  < X < oo 
-00  < X < 00  } 

‘or  a < \ < b, 

lim  P([r(A,n  ,ng)  - g(>v)]/v(A,n1,n2)  < x}  = 4>(x)  , -«•  < x < « . 

n^-> » 

n^oo 

n2^CnJ 

where  c is  an  arbitrary  positive  constant, 

?(*,ni,n2)  = (Tv,-?^)"1  [(^2-A).f(A1,n1)  + (V^)  ?(>y,,n2)] 

and 

v(A,n;L,n2)  = [(7^-Tv)2  fj2(A1,n1)/n1  + Ck-\f  f O/^  VV  . 

Proof.  Let  a double  arrow  denote  convergence  in  distribution  and  let  N(p., 0) 
denote  a normal  random  variable  with  mean  |j  and  standard  deviation  a.  From 
the  conditions  of  the  proposition,  we  may  write 

n^/2[r(X1,n1)  - g^)]  ==£>  N(0,  t ](\)) 

and 

[cn1]1/2[r(>>2,n2)  - g(Ag)]  ==>  N(0,  t]^)) 

as  n^  — •>  00  and  rig  = [cn^] . Since  [cn^]'L/2/(cn1)^2  — > 1 we  have 
(cn1)1/2  [r(*2>n2)  - g(Xg)]  ==>  N(0,  r^Ag)) 

# 

or 


lim  P(n^/2[f(A1>n1)  - gO^)]/^)  < x}  = c(x)  , 


n^-» 00 


and 


lim  Ptng^frCXg^g)  - g(A2)]/ti(Ag)  < x)  = <&(x)  , 


n2-»» 


where  4>  is  the  standard  normal  distribution  function.  Then  f 


k. 


nJ;/2[?(A2,n2)  - g(A2)]  ==>  N(0,  n(>V,)/cl/2)  • 


Since  rCA-^n.^)  and  r(Ag,n2)  are  independent,  we  have,  "by  Theorem  3. 2 
of  [1],  convergence  of  the  vector  process 


nl/2[H\>\)  - g(\)] 


n^2[r(A2,n2)  - g(A2)] 


= ==> 


/ 


N(0,  \ 

\N(0,  nCAgJ/c1/2)  / 


where  the  limit  random  variables  are  independent.  Using  the  Continuous  Mapping 
Theorem,  cf.  [1]  Theorem  5.1,  we  have 


n^2(A2-A1)  1 [(Ag-A)  r(A1,n1)  - (Ag-A)  g(A1)  + (A-A.^  r(yn2)  + (A-A1)  g(Ag)] 
===>  N(0,  [(Ag-A)2  ^(A^  + (A-Ax)2  n2(A2)/c]1/2/(A2-A1)) 


or 


n^2[f(A,n  ,n  ) - g(A)] 

5—5 ±-5^-5 TP? 7 ===>  N(0,1)  . 

[(Ag-A)2  n2(A  ) + (A-A  f ti  (7v,)/c]1/2  (Ag-A  ) 1 


Since  [cn^/c^ ■>  1,  ^(A^n^  ■>  and  n(A2>n2)  > ^(Ag)  in 

probability,  this  implies 


[r^n^ng)  - g(A)] 

[(Ag-A)2  ^2(A1,n1)/n1  + (A-A^2  ^2(A2,n2)/[cn1]]1/2  (Ag-A^'1 


===>  N(0,1) 


or 


5. 


r^n^rig)  - g(A) 
v(X,n~,n2) 


===>  N(0,l) 


which  is  the  desired  result. 


The  above  proposition  allows  one  to  obtain  a (lOO)(l-r)$  confidence 
interval  for  g(A)  as  follows.  For  n^  and  n^  sufficiently  large, 

P(-4>  L(1  - |)  < [r(X,n1,n2)  - g(?\)]/v(A,n1,n2)  < «5_1(l  - |)}  = 1 - Y • 

This  may  be  rewritten  as 

PfrCAjn^ng)  - vC^n^ng)  <b  1(l-  £)  < g(A)  < r(A,n1,n2)  + y(A,n1,n2)<I>  1(l-  £)} 

= i - r, 

giving  the  desired  confidence  interval. 

Note  that  the  confidence  interval  obtained  at 

r(A1,n1)  + n"1/2  r^A^n^  <d"1(1 

which  is  exactly  the  same  confidence  interval  which  could  be  obtained  based  on 
the  n^  observations  at  A = A^ . Similarly,  the  confidence  interval  obtained 
at  A = Ag  is 

r(A2,n2)  + n”1^2  ft^ng)  0_1(l  - |)  . 


A = \ 


reduces  to 


2>  ' 
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Letting  i(A,n^,ng)  denote  the  length  of  the  confidence  interval  obtained 
at  A,  we  see  that 


where  i(A^,n^)  and  are  the  lengths  obtained  at  A^  and  Ag 


respectively.  Hence  i(A,n^,ng)  is  the  square  root  of  a quadratic  function 


3.  Applications  Using  the  Regenerative  Method 

In  this  section,  we  give  examples  illustrating  the  use  of  Proposition  1 
for  a specific  class  of  simulations,  namely,  those  simulations  where  the 
regenerative  method  is  applicable.  (References  [2],  [3],  [4],  [5],  and  [7] 
may  be  consulted  for  more  details  of  the  regenerative  method.)  A basic 
problem  in  statistical  inference  in  simulations  is  to  estimate  the  constant 
g(A)  = E(f(x(A)))  where  f is  a general  real- valued  function,  A is  an  input 
parameter,  and  X(A)  is  the  stationary  random  vector  associated  with 
(X(s>^)'.s  > 0),  the  process  being  simulated.  In  the  regenerative  method,  we 
observe  the  process  (x(s,A):s  > 0}  in  random  cycles  of  lengths  a^(A), 
q^(a),  ...  , C*  (A)  and  record  in  each  cycle  the  values  Y-^A),  Yg(A), . . . ,Y^(A) 
where  Y.(A)  is  the  area  under  the  curve  f(X(s,x))  in  the  i^  cycle. 


a2 

which  passes  through  (A^,  l (A^,n^))  and 


1 


ro|-s 
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The  crucial  conditions  required  for  the  regenerative  method  to  be  used  are 
that  the  n pairs  ((y^(A),  0^(a)),  i = 1,  2,  . . . , n}  are  independent  and 
identically  distributed  and  that  g(A)  = E{ Y^A^/Eta^A)} . Now  define  the 
column  vector  U^A)  = (Y^A),  (^(A))  and  let 


E(A)  = 


ou( A)  o^A)  \ 


a 12(X)  °ZZM  / 


denote  the  covariance  matrix  for  U^(A).  Denote  the  sample  mean  by 


U(A,n)  = 


and  the  sample  covariance  by 


Y(A,n) 
a(A,n)  / 


Z U4(A) 
n 1=1  1 


s(A,n)  = 


sn(A,n)  s^CAjn) 


s21(A,n)  s22(A,n) 


n 


= Z [U  (A)-U(A,n)][U  (A)-U(A,n)]  ' 

n-x  i=1  l l 


where  the  prime  denotes  transpose.  Next,  define  point  estimates  for  g(A) 
as  follows: 

( l)  Classical  estimator 


.<*,«)  - ^ ; 

a(A,n) 


(2)  Beale  estimator 

* , . y(\  n't  f1  + sip(^>n)/n  Y(A,n)  a(A,n)] 

rb(A,n)  = . i± 

a(A,n)  [1  + s22(A,n)/n  a (A,n)l 
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(3)  Tin  estimator 
?.  (A, n)  = 

a(A,n) 

(U)  Jackknife  estimator 

r (A,n)  = ^ Z 0,(A,n) 
J n . , l 


1 + 


s^Ajn) 


Y(A,n)  a(A,n)  a (A,n) 


s22(?v,n)  ! 

12 ' n 


where  0,(A,n)  = n[Y(A,n)/a(A,n)]  - (n-l)  [ Z Y (A)/  Z a.(A)]  . 
1 J j/i  J 


Finally,  define 


nc(A,n)  = [sL1(A,n)  - 2r  (A,n)  s^n)  + r^(A,n)  s22(A,n)]/a(A,n) 

and 

T)T(A,n)  = { Z [0.  (A,n)  - r (A,n)]2/(n-l))1^2  . 

J i=i  1 J 


Now  let 

Zi(A)  = Y.(A)  - g(A)  a. (A) 

p 

and  note  that  E{Z^(A)}  = 0 and  define  a (A)  = var{Z^(A)}.  Since  the  vectors 
(U^( A) , i > i}  are  i.i.d.  it  follows  that  (Z^(A),  i > l)  are  i.i.d.  By  the 
central  limit  theorem  for  partial  sums  of  i.i.d.  random  variables,  it  follows 
that 


lira  P(  Z Z.  (AVn1/2  a(A)  < x}  = $(x)  , -«  < x < » , 

n * i=l  1 


which  may  be  rewritten 


9. 


lira  P{n1^2[r  (X,n)  - g(A)]  a(A,n )/a(A)  < x]  = <t(x)  , -~  < x < » . 

n -*« 


Since  a(X,n)  — > E{G^(X)}  a.e.,  it  follovs  that 

lim  Ptn^^tr  (X,n)  - g(A)]/t|(X)  < x}  = G>(x)  , < x < °o  , 

n C 

where  q(A)  = cj(?0/E[q^(X)}  . Now  it  may  he  shown  that 
nV2[£  (x,n)  - r (X,n)]  ■>  0 a.e., 

C D 

n1/2 [ r (X,n)  - r (X, n)]  ->  0 a.e., 

c z 

n1/2 [ ( A , n ) - rT(X,n)]  > 0 a.e.. 


and  that  tj  (X,n)  — ■>  t|(A)  and  f|T(X,n)  — ■>  q(X)  is  probability  as  n — ■>  ». 

C J 

Hence,  the  conditions  of  Proposition  1 are  satisfied  with  any  of  the  following 
substitutions : 


VW-  or  vw 


substituted  for  k = 1,  2 


°r  V\,nk^  substituted  for  k = 1,  2. 


To  illustrate  the  application  of  Proposition  1,  consider  a simulation  of 
the  customer  waiting  time  process  (W^,  n > 1)  in  an  M/M/l  queue.  Suppose 
we  wish  to  study  the  sensitivity  of  the  mean  stationary  waiting  time 


10. 


g(A)  = E{W(A)}  to  the  arrival  rate  A over  the  range  5 < A < 5,  with 

the  service  rate  ^ = 10.  In  what  follows,  we  shall  illustrate  the  proposition 

for  this  simulation  using  the  "classical"  estimators  r^  and  fjc,  though  we 

could  have  chosen  any  of  the  estimators  given  above . 

In  the  queueing  simulation,  we  say  that  the  i busy  cycle  is  initiated 

vith  the  arrival  of  the  i customer  to  find  an  empty  queue.  Suppose  that 

simulation  runs  consisting  of  nq  = n2  = ^0*000  busy  cycles  are  made  at 

parameter  settings  A = A^  = 3 and  A = A^  = 5.  Let  a^A)  denote  the  number 

th 

of  customers  served  in  the  i busy  cycle,  and  let  Y^(A)  be  the  sum  of  the 
waiting  times  for  those  customers.  Then  it  may  be  shown,  cf.  [2],  that 
((Yi(A),  a^(A),  i > 1)  are  independent  and  identically  distributed,  and 
E{W(A))  = E(Y1(A))/Eta1(A)}.  Hence,  it  is  appropriate  to  apply  the  regenerative 
method  as  discussed  above.  In  particular,  we  can  compute,  for  k = 1,  2, 

Sll^7Vnk^  S12^\,nk^  and  S22^\'\^  and  from  these 
we  can  compute  r (A^o^)  and 

Suppose  that;  as  the  result  of  these  computations, 

9 (A.,n.)  = .042 
c 1 1 

rc(A2,n2)  = .098 

f)c(A1,n1)  = .172 
( Ag , ng ) = -^30 

(Ag-Ai)"1  [(Ag-A)  rc(A1,n1)  + (A-A1)  rc(Ag,n2)] 

.028a  - .042 


Ihen 

rU^n-^ng)  = 


as  defined  above. 


and 
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$(A,tyn2)  = [(X2-X)2  ^(A^)^  + (A-A^2  ^(A^)/^]1/2/^^) 

= [(5-A)2  (.0296)  + (A-3)2  (.l85)]1/2/200  . 

Now  assume  that  E(W(A))  is  approximately  linear  for  3 < A < 5 and 
suppose  that  we  desire  a 98$  confidence  interval  for  E(W(A))  for  each  A 
such  that  3 < A < 5*  From  Proposition  1,  such  a confidence  interval  is  given 
by 

fCAjn^ng)  + v(X,n1,n2)  $_1(l  - |)  , 

where  r(A,n^,n2)  and  v(A,n^,n2)  are  computed  as  above  and  y = .02.  These 
confidence  intervals  are  illustrated  in  Figure  1.  Thus,  for  example,  a 98$ 
confidence  interval  for  E{W(U)}  is  .070  + .005^. 

It  is  of  interest  to  compare  the  confidence  intervals  shown  in  Figure  1 
to  a confidence  band  obtained  for  the  function  E(W(A)}  using  the  method  of 
[2].  Figure  2 shows  such  a comparison.  A 96 ‘jo  confidence  band  for  E(W(A)} 
over  3 < A < 5 consists  of  two  straight  lines  intersecting  the  98$  confidence 
intervals  at  A = 3 and  A = 5 and  everywhere  containing  the  98$  confidence 
intervals  at  3 < A < 5*  Urns,  if  one  desires  confidence  bands  rather  than 
confidence  intervals,  it  is  apparent  that  such  confidence  bands  are  obtained 
at  a cost  of  lower  levels  of  confidence  in  addition  to  wider  intervals . 

Suppose  that  after  running  simulations  at  A = 3 and  A = we  are 
interested  in  estimating  E(W(A))  for  2 < A < 5*  If  we  assume  that 
E{W(A)}  is  approximately  linear  over  2 < A < 5>  then  the  above  expression 
for  the  confidence  intervals  remains  valid.  See  Figure  3 f°r  an  illustration 
of  the  confidence  intervals  so  obtained. 

Acknowledgement . It  is  a pleasure  to  acknowledge  useful  discussions  with 
Donald  L.  Iglehart  concerning  the  research  reported  here. 


Mean  Waiting  Time,  E(W(A)} 


12. 


Function  of  the  Arrival  Rate  A Over  3 < A < 5« 
(m/m/1  Queue  with  Service  Rate  ^ = 10. ) 


13. 


Figure  2.  Comparison  of  Confidence  Intervals  vs.  Confidence  Band  for 
the  Mean  Waiting  Time,  as  a Function  of  the  Arrival  Rate  X 
Over  3 < X < 5* 


Mean  Waiting  Time,  E{w(A)) 


14. 


Fig.  3.  98$  Confidence  Intervals  for  the  Mean  Waiting  Time  as  a 

Function  of  the  Arrival  Rate  A Over  2 < A < 5* 
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In  this  paper,  we  continue  the  investigation  of  response  surface 
methodology  in  simulation,  begun  in  references  [2],  [3],  and  [6].  Suppose 
that  g is  a linear  output  function  of  an  input  parameter  A,  over  some 
range  a < A < b,  and  that  observations  are  taken  at  two  parameter  settings 
Ai  and  It  is  shown  that  if  certain  conditions  are  satisfied,  then  a 

confidence  interval  may  be  obtained  for  g(A)  at  any  point  a < A < b 
without  making  any  additional  observations . An  Illustration  is  provided 
which  makes  use  of  the  regenerative  approach  for  analyzing  simulations. 
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